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Abstract. We present a new spectral-method-based algorithm for finding ap- 
parent horizons in three-dimensional space-like hypersurfaces without symmetries. 
While there are already a wide variety of algorithms for finding apparent horizons, 
our new algorithm does not suffer from the same weakness as previous spectral 
apparent horizon finders: namely the monopolar coefficient (£ = in terms of the 
spherical harmonics decomposition) needed to be determined by a root-finding 
procedure. Hence, this leads to a much faster and more robust spectral appar- 
ent horizon finder. The finder is tested with the Kerr-Schild and Brill-Lindquist 
data. Our finder is accurate and is as efficient as the currently fastest methods 
developed recently by Schnetter (2003 Class. Quantum Grav. 20, 4719) and 
Thornburg (2004 Class. Quantum Grav. 21, 743). At typical resolutions it takes 
only 0.5 second to find the apparent horizon of a Kerr-Schild black hole with 
a = 0.9M to the accuracy ~ 10~^ for the fractional error in the horizon's location 
on a 2 GHz processor. 



PACS numbers: 04.25.Dm, 04.70.Bw, 02.70.Hm 

1. Introduction 

Apparent horizons play an important role in numerical relativity for spacetimes 
containing black hole(s). Being defined locally in time (see section [2]), the apparent 
horizon(s) can readily be computed from the data on each hypersurface during a 
numerical evolution in (3 -I- 1) numerical relativity. In contrast, the event horizon is 
a global property and can be determined approximately only when the spacetime has 
essentially settled down to a stationary state. Once the spacetime has settled down, 
the event horizon can be found at all previous times by integrating null geodesies 
backwards in time (e.g., [Ij [H [S] |4] ) . As it (if exists) must be inside an event horizon 
[5], the apparent horizon is an important tool to track the location and movement 
of black hole(s) in a numerically generated spacetime. Furthermore, the surface of 
an apparent horizon also provides a natural boundary within which the spacetime 
region can be excised from the computational domain in order to handle the physical 
singularities inside a black hole O [71 [8] (also see, e.g., [9l [TOl [TTl [12] for black hole 
simulations without excision). In the new concept of a 'dynamical horizon' (see 
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[131 [T3] for reviews), apparent horizons are essentially the cross sections of the (three- 
dimensional space-like) dynamical horizon on the hypersurfaces. It has recently been 
shown in this context that the areas of the apparent horizons satisfy a causal evolution 
equation and give a positive bulk viscosity in a viscous fluid analogy [15|. This is in 
contrast to the event horizon which yields a noncausal evolution and a negative bulk 
viscosity. 

A wide variety of algorithms for finding apparent horizons have been proposed 
in the past decade. We refer the reader to the review article by Thornburg [16] (and 
references therein) for details. In this paper, we present a new apparent horizon finder 
which is based on spectral methods. While the spectral-method-based algorithm for 
finding apparent horizons was first proposed by Nakamura et al [17] more than 20 years 
ago, our new approach does not suffer from the same weakness as in the Nakamura et 
al algorithm: namely the £ — coefficient of the spherical harmonics decomposition of 
the apparent horizon's surface needed to be determined by a root-finding procedure. 
Hence, our algorithm leads to a more robust and efficient spectral apparent horizon 
finder. We have tested our finder with analytic solutions for single and two black-hole 
spacetimes. Our finder is as efficient as the currently fastest algorithms developed by 
Schnetter [18] and Thornburg [T9]. 

This paper is organized as follows. In section [2| we present the notations and 
various definitions. In section [3| we briefly review the Nakamura et al. algorithm; we 
describe our spectral algorithm and the numerical procedure in section |4| Section [5] 
presents tests with analytic solutions to assess the accuracy, robustness, and efficiency 
of our finder. Finally, we summarize our results in section [6l Latin (Greek) indices go 
from 1 to 3 (0 to 3). 



2. Notations and definitions 



Given a space- like hypersurface E with future-pointing unit normal n'^, the 3-metric 
7^1^ induced by the spacetime metric g^j/ onto E is 

Ifj^i- ■= 9t^u + n^n,,. (1) 

Let S* be a closed smooth (two-dimensional) surface embedded in E, and let be 

the outward-pointing unit normal of S*, which is spacelike and also normal to (i-e., 
Sps'' = 1 and s^n^ = 0). The 3-metric now induces a 2-metric on S: 

Ifiv - SpSi.. (2) 



Let k'^ be the tangents of the outgoing future-pointing null geodesic whose projection 
on E is orthogonal to S. We have (up to an overall factor) 

fc^ = s^-|-n'", (3) 

on the 2-surface S. 

The expansion of the outgoing null geodesies is 

e = v^fc^ (4) 

where is the covariant derivative associated with g^^. In terms of three-dimensional 
quantities, on the 2-surface S, the expansion can be written as (see, e.g., [20] ) 

e = As* - K + s's^K,^, (5) 

where Di is the covariant derivative associated with 7^ , Kij is the extrinsic curvature 
of E and K is the trace of Kij . The expansion can also be written as 

e = m'^ {D,sj - K,,) . (6) 
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The 2-surface S is called a marginally trapped surface if 8 = everywhere on S. We 
shall call here the outermost of such surfaces (which is a marginally outer trapped 
surface - MOTS) the apparent horizon. 

To parameterize the apparent horizon, we assume that the topology of S* is a 
2-sphere, and S is star-shaped around the coordinate origin r — 0, which means that 
for every point M inside S, the straight line connecting the origin to M is entirely 
inside S [21]. The position of the apparent horizon can then be represented as 

F{r,9,ip):=r-h{9,ip) = 0, (7) 

where (r, 9, tp) are the standard spherical coordinates. The function h measures 
the coordinate distance to the horizon's surface in the direction {9,Lp). With this 
parametrization, the unit normal s* is given by 

D^F D'F 



{r'D,FD,Ff'^' \DFy 
where := Dj. The expansion (equation ([5])) becomes 



(8) 



where the condition m^^ Sj = has been used. 



3. The Nakamura et al algorithm 



In this section, we give a brief review of the algorithm adopted by Nakamura et al |17| 
for finding apparent horizon based on spectral methods. They expand h in spherical 
harmonics: 

h{9,^) - £ ^ a,^Yr{9,^). (10) 

1=0 m=-l 

They rewrite the apparent horizon equation 8 = as|]j] 

Ag^h ^ pQ + Ag^h, (11) 

where Ag^ is the flat Laplacian operator on a 2-sphere defined by 

Ag^h := hfie + cot 9hfi + sin"^ 9h^ipip. (12) 

The positive scalar function p is chosen such that the term h^g cancels on the right- 
hand side (RHS). Using the fact that the F/" are an orthogonal set of eigenfunctions 
of Ag^, 

Ag^vr = -l{l+\)Yr. (13) 
we obtain the relation (with d£l = sin i 



£ (^ + 1) a,™ = / Vr (p8 + Ag^h) dVt. (14) 



s 



This equation can be used to solve for the coefficients aim via an iteration procedure. 
However, the value of aoo has to be determined at each iteration step by solving for 
the root of 

/ Fo°* (P0 + Ae^/i) d^l = 0. (15) 
Js 

II We follow the notation of Gundlach 1221. 
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The main disadvantage of the above scheme is that the coefficient ago has to 
be determined separately by equation (fT5|l. As pointed out by Gundlach [22, solving 
equation (jl5p by any iteration method is as computationally expensive as many steps of 
the main iteration loop. Furthermore, equation (jlSp may have multiple roots or none. 
In those cases, each root or each minimum (if there is no root) should be investigated 
separately [23] . This clearly reduces the efficiency of the algorithm significantly. 



4. Our algorithm 



Master equation for apparent horizon 

Our spectral-mcthod-based algorithm uses a similar ansatz (jlip as Nakamura et al 
|17| . The main difference is that we do not need to determine aoo separately. Hence, 
this leads to a more robust and efficient apparent horizon finder based solely on the 
spectral methoc0. 

To begin, we first introduce a flat metric fij on the hypersurface S. The 
components of the flat metric with respect to the spherical coordinates {r,9,ip), and 
the associated natural basis (^, ^), a-re fij = diag(l, r^, sinfl). Let Di be the 
covariant derivative associated with fij. The expansion function Q (equation ([9|)) can 
now be written as 

e = {Y' - s's') [\DF\-^ (DiDjF - A'^^jD^F) - K,j\ , (16) 

where the tensor field A™j is defined by 

We have also used the following relation between the two covariant derivatives Di and 
Di-. 

D^Vj ^D,Vj ~ A"ljV^, (18) 

where is an arbitrary 3- vector on E^. 

Motivated by the recently proposed fully constrained-evolution scheme for 
numerical relativity [24j . we define a conformal factor ^! by 

and also a tensor field /i*-' by 

= ^-4 ^ f^^3^J (20) 

We also expand all tensor fields onto the following spherical basis: 

e — — e- — i— c — ^ ^ (21) 
Qj,' r 86' r sin d dip 

This basis is orthonormal with respect to the fiat metric: = diag(l, 1,1). Here and 
afterwards, we denote the tensor indices associated with this basis with a hat. The 
expansion function now becomes 

e = ■^-^\DF\-^ f'^D{D-.F + {^^-'^l^'^ - s*s^) \DF\-^D{D.F 



1 See 1221 for Gundlach's "fast flow" algorithm which combines the spectral algorithm of Nakamura 
et al and the so-called curvature flow method (see also | 16| for discussion). 
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_ (yi _ s^J^ {\DF\-'A^.D,nF + K^^ . (22) 
Now let us consider the first term on the RHS of this equation: 
^-^\DF\-^fW{D.F = ^-^\DF\-^ (Df^DfF + DgDgF + D^D0F) 

= T-An^ni 9 i^^,ee + coiOhfi + sin"^ 9h,u,u, - 2r) 
^^\DF\r^ ^ I 

where we have set r = h(Q^ ip) for the apparent horizon in the last equality. In 
equation Ij23[) , we have used the following relation for the components of the covariant 
derivative -D-- of a 3-vector in the orthonormal basis {ej}: 

^^-^/a^^-r^i^' (24) 

where ef' := diag[l, 1/r, l/(rsin6')]. The F'^,-. are the connection coefficients of Z?j, 
associated with {ej}. The non- vanishing components are 

yf _ _ yf _ flip _ _ _yV _ ~^ (0K\ 

68 f6 ^ v0 ^' ^ w rtane' ^ ' 

Equations ((22)l and ((23| suggest that, instead of the ansatz (fTTj) as taken by 
Nakamura et aZ, it is more appropriate to rewrite the apparent horizon equation = 
as 

Ae^/i -2h = \Q + Ae^h - 2h, (26) 

where the scalar function A is chosen to be A = ^''*|Z3i^|/i^ such that the combination 
Agj^/i — 2h cancels on the RHS of this equation. Hence, the master equation that we 
solve in our algorithm is 



{^DF\-^A%praF-rK:^\. (27) 



The expansion coefficients are now determined by solving the following equation 
iteratively: 

= ^ ^ ^ ^ Vr (A9 + Ae^fe - 2h) dn. (28) 

This equation applies for all > 0, and hence ago is not treated specially. In [25], 
Shibata developed an apparent horizon finder based essentially on the same form of 
equation p6p . However, he solved the equation using the finite-differencing method 
without pointing out the key advantage that, if solved by the spectral method, 
the coefficient ooo (as determined by our equation (|28[) and his equation (1.3) in 
[25] ) does not needed to be solved by the root-finding procedure. In this work, 
we solve the algorithm for the first time with spectral method. The difference 
between equations ([25]) and IH]) leads to a dramatic improvement in the efficiency 
and robustness of spectral-method-based algorithms for finding apparent horizons. 
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4-2. Numerical procedure 

For given 3-metric and extrinsic curvature Kij on a hypersurface E, equation (|27p 
represents a nonlinear elliptic equation for the function h. We solve this equation 
iteratively by considering the RHS of the equation as a source term for the linear 
operator Ag^p — 2 acting on h. We use a multidomain spectral method to solve the 
elliptic equation [211 [26]. The code is constructed upon the C++ library LORENE 
[27] . and is publicly available. 

The numerical iteration procedure is briefly described here. Assume that the 
data i'^ipK^) are given on E. The conformal factor 5* and the tensor field K^^ are 
then calculated by equations and ([^0]) , respectively. Assume that an initial guess 
for the function h{9,ip) is chosen (equivalently for the spectral coefficients aim)- The 
iteration processes are as follows. 

(i) At the nth iteration step, the function ft,*^"^ is determined by the coefficients 
(with the superscript (n) labels the iteration steps) . The level-set function F and 
the unit normal vector are then obtained from /i^") (see section [2]). 

(ii) The spectral coefficients at the next iteration step are calculated by equation ([^5]) : 

- ... , o / rrs^-^dn, (29) 



^(^+l) + 2 



s 



where S''-"-' represents the RHS of equation ([27]) evaluated from a^"). The new 



function is then obtained from a^^^^'' by equation ([TO 



(iii) The difference between and /i^"-' is calculated. The iteration procedure 

continues until the maximum value of the difference throughout the whole angular 
grid (9i , ipj ) is smaller than some prescribed value eh ■ 



5. Tests 



5.1. Kerr-S child data 

As a first test of the apparent horizon finder, we use a single black hole in Kerr-Schild 
coordinates (see, e.g., jS^) to study its convergence properties and robustness. Let 
M and a denote respectively the mass and spin parameter of the black hole. In the 
standard spherical coordinates (r,9,ip)^ the polar and equatorial coordinate radii of 
the apparent horizon are given by 

rpo = r, Tcq = Vr'^ + a^, (30) 

where r = M + \/ — a^ . The area is given by 

A = 4:TT (r^ + a^) . (31) 

We first test the convergence property of the code with respect to increasing 
number of collocation points from runs with a black hole of M = 1 and a = 0.9. 
The polar radius of the apparent horizon is rpo ~ 1.436 and the equatorial radius is 
Tcq ~ 1.695. The analytic data (Tij^^ij) set on the numerical grid points of the 
computational domain ranging from r = 1 to r = 5, which is covered by three spectral 
domains. The boundary between the first and the second domains is at r = ri2 = 1.5, 
whereas that between the second and the third domains is at r = r23 = 2.5. In each 
domain, we use {Nr, Ng, N^p) collocation points. We also enforce a symmetry with 
respect to the equatorial plane. The initial guess for /i is a sphere at r = 3. 
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Table 1. Convergence test for a Kerr-Schild black hole with M = 1 and 
a = 0.9. Listed are the number of radial collocation points in each domain Nr, the 
fractional errors in the polar (equatorial) coordinate radius Arpo/rpo (Arcq/roq) 
and area A A/ A, the maximum remaining error of the expansion function ASmax 
on the horizon and the run times. We use Ng = (A'^r + l)/2 points in the polar 
direction and = 1 in the azimuthal direction. We set the iteration parameter 



ejl = 


10-10. 


















Nr 






A^'cq/ Tcq 




AA/A 




A0niax 




Time (s) 


13 


2.360 X IQ- 


-5 


8.003 X 10- 


~B 


1.276 X 10" 


-b 


4.240 X 10" 




0.507 


17 


1.693 X 10- 


-6 


9.705 X 10- 


-8 


1.605 X 10- 


-6 


5.333 X 10- 


4 


0.747 


21 


1.580 X 10- 


-7 


2.145 X 10- 


-8 


1.788 X 10- 


-7 


6.574 X 10- 


5 


1.129 


25 


1.692 X 10- 


-8 


4.235 X 10- 


-9 


2.045 X 10- 


-8 


8.033 X 10- 


6 


1.615 


33 


1.067 X 10- 


10 


1.707 X 10- 


10 


4.733 X 10" 


10 


1.559 X 10- 


7 


3.059 


37 


1.649 X 10- 


10 


1.325 X 10- 


10 


2.863 X 10" 


10 


2.383 X 10- 


8 


4.286 


41 


1.590 X 10" 


10 


1.053 X 10" 


10 


2.154 X 10" 


10 


3.790 X 10" 


9 


5.722 



Table [T] shows the results for increasing Nr, with A''^ = (A^r- -l- l)/2 and N^p = 1. 
We choose the iteration parameter eh = 10~^° in this test (see section IT2|) . In the 
table, for each Nr, we list the fractional errors in the polar (equatorial) coordinate 
radius Arpo/rpo (Aroq/roq) and area AA/A, the maximum remaining error in the 
expansion function AO^ax on the horizon's surface, and the run time^U- The error in 
the area is defined by AA/A := |(^ana ^ ^num)Mana|, where the analytic result Aana 
is given by equation pip and the numerical result is calculated by the integral 

y/^h^ sin edOdip, (32) 

's 

with q being the determinant of the 2-metric on the apparent-horizon's surface 
(expanded onto the basis {e^}). Explicitly, in terms of a general 3-metric 7j-, Annm is 
given by 



2tt ^tt 



^num = / / [{iffh g + 2j-ghh,g + jggh'^) [iffh^ + sin + •^^^h 

Jo Jo 



^ sin^ ( 



,ip~^Jfiphh^gsin6-\-'fg^h sin^j 



-,1/2 

dOdip. (33) 



In figure[T]we plot AA/A against Nr to explicitly show the convergence behaviour 
of the finder for three different choices of eh- It can be seen that the error AA/A 
converges exponentially towards zero with the number of points, as expected for 
spectral methods, until the accuracy is limited by the choice of eh- Furthermore, we 
also see that the number of iterations to a given error level eh is essentially independent 
of the value of £max used in equation pO]) . This agrees with the conclusions obtained 
from the original Nakamura et aPs algorithm (or its modifications) as investigated by 
Kemball and Bishop ^23j . 

Next we test the robustness of our finder by performing runs with different initial 
guesses for h- We use the same black hole as above (M = l,a = 0.9), but with a 
larger computational domain ranging from r = 1 to r = 10. The boundaries between 
the different spectral domains are ri2 — 2.5 and r23 = 5.5. In general, we set up an 

The run times correspond to the CPU time the code took to locate the apparent horizon on a 2 
GHz Intel Core Duo processor. The best-fitted curve suggests that the scaling of the run time is close 
to Ng, which comes from the computation of discrete Legendre transforms. 
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Figure 1. Convergence towards zero of the fractional error in the area AA/A 
with the number of collocation points for three different choices of the iteration 
parameter e^. 

Table 2. Robustness test for the same black hole used in table [T] The polar 
coordinate radius of the apparent horizon is rpo ~ 1.436 and the equatorial 
coordinate radius is rcq ~ 1.695. The initial guess for the 2-surface h is given 
by the surface of an ellipsoid with axes (a, b, c) defined in equation I I34I I. We use 
collocation points A^r = 25 and Ng = 13 for all the five cases considered, but 
= 1 (4) for cases A and B (C-E). 



Case 


(a, b, c) 


AA/A {10-") 


A 


(8,8,8) 


4.04871 


B 


(1.2, 1.2, 1.2) 


4.04905 


C 


(4,6,8) 


4.04875 


D 


(2,3,1.2) 


4.04910 


E 


(1.2,1.5,2) 


4.04909 



initial guess for h to be the surface of an ellipsoid given by 

? 2 2 

- + I2 + - = 1' (34) 
0^ c"^ 

where {x,y,z) are the Cartesian coordinates relating to the spherical coordinates 
(r, 0, if) in the standard way. The constants (a, 6, c) are freely chosen. The initial 
guess used for the results listed in table [1] corresponds to a = b = c = 3. Table [2] 
contains the results for five different initial guesses. Case A corresponds to a sphere 
with a coordinate radius which is about five times away from the horizon's surface. 
On the other hand, the initial surface for case B is a sphere located entirely inside the 
apparent horizon. The initial guess for case C is an ellipsoid enclosing the horizon. 
Finally, cases D and E represent initial surfaces which cross the horizon. The results 
show that our finder can locate the apparent horizon to the same accuracy with all 
the five (quite generic) choices of (a, 6, c). 

One of the main requirements of an apparent horizon finder is speed. This is in 
particular an important issue if the finder has to run frequently during a simulation. 
In order to compare the speed of our finder with some other commonly used methods, 
we take the data given by Schnetter [TH]. 
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Table 3. Schwarzschild black hole offset from the coordinate origin. The hole 
is located at the Cartesian coordinates {d/\/2,d/\/2,0). Listed are the offset d, 
the fractional error in the area AA/A, and the maximum remaining error of the 
expansion function on the horizon AQmax- 



d 


AA/A 




0.1 


9 X 10-" 


3 X 10- 


-a 


0.2 


3 X lO-** 


5 X 10" 


-3 


0.3 


2 X lO"** 


2 X 10" 


-2 


0.4 


1 X 10-"' 


5 X 10" 


-2 



In table 5 of T8] Schnetter compared the run times to locate the apparent horizon 
of a Kerr-Schild black hole with M = 1 and a = 0.6 for his elliptic method and two 
other methods, namely the fast-flow 22J and minimization [29] algorithms. The fastest 
case (0.5 s on a 1.2 GHz processor) was obtained by his elliptic method with the initial 
guess being a sphere at r = 2. The error in the area is AA/A = 9 x 10"'^ (according 
to table 4 of 18J). For comparison, the fast-flow and minimization algorithms took 
more than 10 s and 90 s respectively in the test [18_. 

We have performed tests with the same black hole and initial guess, and found that 
our finder took 0.129 s (on our 2 GHz processor) to locate the horizon to the accuracy 
AA/A = 2 X 10"^ using the resolution {Nr,Ng,N^) = (7,5,1) with eh = 10"^. We 
have also used Thornburg's finder AHFinderDirect [1£^ (which is implemented 
within the CACTUS computational toolkit [30]) to perform the same test using 
Cartesian grid resolutions {N^,Ny,N^) = (31,31,19) with Aa; = Ay = Az = 0.2 
in bitant symmetry. We found that his finder took 1.004 s (on our 2 GHz processor) 
to locate the apparent horizon to the accuracy AA/A = 3 x lO"**. 

We note that the above test does not represent a direct comparison between 
the different algorithms because of the different grid structures (Cartesian versus 
spherical coordinates), code implementations, memory usage and computer systems. 
Nevertheless, we can conclude that for this particular test, to obtain about the same 
accuracy level, our spectral-method-based finder is as efficient as the finders developed 
by Schnetter [1S| and Thornburg [TO] . 



5.2. Brill- Lindquist data 

In this part, we test our finder using the Brill-Lindquist data [31]. This is a classic test 
involving multiple black holes used in numerical relativity. The 3-metric is conformally 
flat, 7y = (f>'^fij, and is time symmetric (i.e. Kij — 0). For two black holes, (f> is given 

by 

2\r-ri\ 2|r-r2| 

where Mi (i = 1, 2) is the mass of the ith black hole and are the coordinate positions 
of the holes. 

We first begin with a single black hole (M2 — 0), in which case the problem 
is equivalent to a Schwarzschild black hole in isotropic coordinates offset from the 
coordinate origin. The apparent horizon is a coordinate sphere of radius Mi/2 with 
respect to the center of the hole. The area of the horizon is A = IGirMf. We set 
Ml = 1 and the coordinate position of the hole at ri = {xi,yi,zi) = {d/y/2,d/y/2,Qi). 
We have varied d in order to verify that our finder also works when the center of 
the spherical harmonics is offset from the center of the horizon. Table [3] lists the 
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results for four different values of d. The initial guesses are always a = 5 = c = 1 in 
equation ([M)) . We use three spectral domains to cover the spatial shce up to r — 1.5, 
with collocation points {N^, Ng, Nif,) = (33,17,16) in each domain. The boundaries 
between the domains are ri2 = 0.5 and = 0.8. Similar to [53], we see that 
the accuracy drops quite significantly for very distorted surfaces with respect to the 
coordinate origin. In particular, the error in the area AA/A increases by almost two 
orders of magnitude when d increases from 0.3 to 0.4; this error could be reduced 
using higher grid resolutior0. Nevertheless, it is worth to point out that the original 
Nakamura et al spectral algorithm [17j would not produce any results for d = 0.3 and 
0.4 because equation ([15]) has no roots [23]. We also see that the results are essentially 
independent of the direction of the offset. 

Next we turn to Brill-Lindquist data for two black holes of equal mass. In 
particular, we take Mi = M2 = 1 in the test. The data forms a one-parameter 
family parameterized by the coordinate separation d between the holes. When they 
are far apart, each hole has an individual apparent horizon. For small separation, there 
is a single common apparent horizon enclosing both holes. Determining the critical 
separation at which the common horizon appears in this two black-hole spacetime is a 
standard test problem for apparent horizon finders. The critical separation obtained 
originally by Brill and Lindquist is dc — 1.56 |31j . while more recent results suggest 
that 4 « 1.53 (e.g., [H [H [31 [SS])- In particular, we note that Thornburg ^ 
and Shoemaker et al [33j report very close results at dc — 1.532 and dc = 1.535 
respectively. Nevertheless, Thornburg reports A — 196.407 for the area of the critical 
apparent horizon, which is quite different from the value A — 184.16 obtained by 
Shoemaker el al.. 

Here we test our finder by trying to find a common horizon at the critical 
separations, as reported by Thornburg 19J and Shoemaker et al. [33j. The black 
holes are on the z-axis, with their centers at z = zLd/2. In the test, we use four 
spectral domains to cover the spatial slice up to r = 2. The boundaries between 
the domains are ri2 = 0.5, r23 = 1 and r^^ — 1.5. The initial guesses are 
a = 6 = c = 2in equation ([31)1 . We use {Nr, Ne, N^) = (41,31, 1) in each domain 
and the iteration parameter eh ~ 10^^. Our finder reports a common horizon at 
d — 1.532 (Thornburg's critical value) with the area of that horizon determined to be 
A = 196.417, which agrees to Thornburg's value to 0.005%. The maximum remaining 
error of the expansion function on the horizon is AG^ax ~ 7 x 10^**. The finder took 
59.2 s to locate the horizon. We note that increasing eh to 10~^ would reduce the run 
time to 23.3 s, without changing the three significant figures of A. On the other hand, 
for the same grid setting and parameters, our finder does not find a common horizon 
at the critical value d = 1.535 reported by Shoemaker et al.. In general, for d > 1.532, 
we find two disjoint apparent horizons surrounding fi and r2 by setting the coordinate 
origin for the apparent horizon finder (see section [2]) separately at around the points 
fi and r2 ■ In figure [5] we show the position of the common apparent horizon on the 
x-z plane for the case d = 1.532. The results obtained by four different values of Ne 
(with Nr — 41 and A^,^ = 1 fixed) are plotted together to show the convergence of the 
horizon. 

* The error A.A/A drops down to 2x 10~^ for the case d = 0.4 using collocation points {Nr, Ng, Ncp) = 
(33,25,24). 
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Figure 2. Position of the common apparent horizon on the x-z plane for Brill- 
Lindquist data with d = 1.532. The black holes are centered at 2 = ±d/2 along 
the z-axis. The results obtained by four different values of Ng are shown (with 
Nr = 41 and A^^ = 1 fixed). 



6. Conclusions 

In this paper we liave presented a new apparent liorizon finder wliicli is based on 
spectral methods. Our proposed algorithm does not need to treat the i = coefficient 
of the spherical harmonics decomposition separately as required in previous spectral 
apparent horizon finders [T71 [53]. Hence, this leads to a faster and more robust 
finder based solely on spectral methods. We have made a performance comparisons 
with other apparent horizon finders using the Kerr-Schild data. Our finder is much 
faster (by orders of magnitude) than other commonly used methods (e.g., the fast-flow 
and minimization algorithms). It is also as efficient as the currently fastest methods 
developed recently by Schnetter [TB] and Thornburg [T^]. We have also shown that 
our finder is capable of locating the horizon of a shifted Schwarzschild black hole with 
a large offset from the coordinate origin. This would not be possible by using the 
original Nakamura et al spectral algorithm p/7j because equation has no roots if 
the offset is too large. We have also tested our finder for a two black-hole spacetime 
using the Brill-Lindquist data. In particular, we have verified previous results on the 
critical separation at which a common horizon appears in this spacetime. 

Our apparent horizon finder is implemented within the C-l — h library LORENE for 
numerical relativity [27j . and is freely available. The finder should be easily adopted in 
spectral- method-based evolution codes [MUMHSSISSIj particularly to those using shell- 
like domains in spherical coordinates. Comparing to other freely available apparent 
horizon finders which are based on the finite-differencing method (AHFinder [32] and 
AHFinderDirect [19]), our finder also represents another option available to finite- 
differencing evolution codes, with some interpolation to be implemented between the 
finite-difference and spectral grids as, for example, in [3T. 

Note added: After we have submitted the paper, Tsokaros and Uryu informed us that 
they had recently developed an apparent horizon finder based on the same formulation 
presented in this paper |38| . Their work and ours were done independently. 
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